SPECTRAL METHODS FOR PARAMETERIZED MATRIX 

EQUATIONS 

PAUL G. CONSTANTINE* , DAVID F. GLEICH 1 " , AND GIANLUCA IACCARINO* 

Abstract. We apply polynomial approximation methods — known in the numerical PDEs 
context as spectral methods — to approximate the vector-valued function that satisfies a linear 
system of equations where the matrix and the right hand side depend on a parameter. We derive 
both an interpolatory pseudospectral method and a residual-minimizing Galerkin method, and we 
show how each can be interpreted as solving a truncated infinite system of equations; the difference 
between the two methods lies in where the truncation occurs. Using classical theory, we derive 
asymptotic error estimates related to the region of analyticity of the solution, and we present a 
practical residual error estimate. We verify the results with two numerical examples. 
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1. Introduction. We consider a system of linear equations where the elements 
of the matrix of coefficients and right hand side depend analytically on a parame- 
ter. Such systems often arise as an intermediate step within computational methods 
for engineering models which depend on one or more parameters. A large class of 
models employ such parameters to represent uncertainty in the input quantities; ex- 
amples include PDEs with random inputs [21 031 [29] , image deblurring models [8] , 
and noisy inverse problems [7] . Other examples of parameterized linear systems oc- 
cur in electronic circuit design p2j , applications of PageRank [3 [9] , and dynamical 
systems [10] . Additionally, we note a recent rational interpolation scheme proposed 
by Wang et al. [27] where each evaluation of the interpolant involves a constrained 
least-squares problem that depends on the point of evaluation. Parameterized linear 
operators have been analyzed in their own right in the context of perturbation theory; 
the standard reference for this work is Kato [21]. 

In our case, we are interested in approximating the vector-valued function that 
satisfies the parameterized matrix equation. We will analyze the use of polynomial 
approximation methods, which have evolved under the heading "spectral methods" 
in the context of numerical methods for PDEs [3] [6] [20] . In their most basic form, 
these methods are characterized by a global approximation of the function of in- 
terest by a finite series of orthogonal (algebraic or trigonometric) polynomials. For 
smooth functions, these methods converge geometrically, which is the primary reason 
for their popularity. The use of spectral methods for parameterized equations is not 
unprecedented. In fact, the authors were motivated primarily by the so-called poly- 
nomial chaos methods [HI [29] and related work [2j [1] [28] in the burgeoning field of 
uncertainty quantification. There has been some work in the linear algebra commu- 
nity analyzing the fully discrete problems that arise in this context [TH [231 [TT] , but 
we know of no existing work addressing the more general problem of parameterized 
matrix equations. 

There is an ongoing debate in spectral methods communities surrounding the rela- 
tive advantages of Galerkin methods versus pseudospectral methods. In the case of pa- 
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Table 1.1 

We attempt to use a consistent and clear notation throughout the paper. This table details the 
notational conventions, which we use unless otherwise noted. Also, all indices begin at 0. 



Notation 


Meaning 


A(s) 


a square matrix- valued function of a parameter s 


b(s) 


a vector- valued function of the parameter s 


A 


a constant matrix 


b 


a constant vector 


(■) 


the integral with respect to a given weight function 


(•)„ 


the integral (•) approximated by an n-point Gauss quadrature rule 


[M] rxr 


the first r x r principal minor of a matrix M 



rameterized matrix equations, the interpolatory pseudospectral methods only require 
the solution of the parameterized model evaluated at a discrete set of points, which 
makes parallel implementation straightforward. In contrast, the Galerkin method re- 
quires the solution of a coupled linear system whose dimension is many times larger 
than the original parameterized set of equations. We offer insight into this contest by 
establishing a fair ground for rigorous comparison and deriving a concrete relationship 
between the two methods. 

In this paper, we will first describe the parameterized matrix equation and char- 
acterize its solution in section [2] We then derive a spectral Galerkin method and a 
pseudospectral method for approximating the solution to the parameterized matrix 
equation in section[3l In section[4j we analyze the relationship between these methods 
using the symmetric, tridiagonal Jacobi matrices - techniques which are reminiscent 
of the analysis of Gauss quadrature by Golub and Meurant ;T7] and Gautschi [T3] . We 
derive error estimates for the methods that relate the geometric rate of convergence 
to the size of the region of analyticity of the solution in section 03 and we conclude 
with simple numerical examples in section [6l See table Q] for a list of notational con- 
ventions, and note that all index sets begin at to remain consistent with the ordering 
of a set of polynomials by their largest degree. 

2. Parameterized Matrix Equations. In this section, we define the specific 
problem we will study and characterize its solution. We consider problems that depend 
on a single parameter s that takes values in the finite interval [—1,1]. Assume that 
the interval [—1, 1] is equipped with a positive scalar weight function w[s) such that 
all moments exist, i.e. 

(s k ) = J s k w{s) ds < oo, A; = 1,2,..., (2.1) 

and the integral of w(s) is equal to 1. We will use the bracket notation to denote an 
integral against the given weight function. In a stochastic context, one may interpret 
this as an expectation operator where w(s) is the density function of the random 
variable s. 

Let the K^-valued function x(s) satisfy the linear system of equations 

A(s)x(s) = b(s), se[-l,l] (2.2) 

for a given R Nx ^-valued function A(s) and R^-valued function b(s). We assume that 
both A(s) and b(s) are analytic in a region containing [—1, 1], which implies that they 
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have a convergent power series 

A(s) = A + Ais + A 2 .s 2 + --- , b(s) = b + b lS + b 2 s 2 H , (2.3) 

for some constant matrices Aj and constant vectors bj . Additionally, we assume that 
A(s) is bounded away from singularity for all s E [—1,1]. This implies that we can 
write x(s) = A^ 1 (s)b(s). 

The elements of the solution x(s) can also be written using Cramer's rule [2"31 
Chapter 6] as a ratio of determinants. 

*<w = rar •=» <"> 

where Ai(s) is the parameterized matrix formed by replacing the ith column of A(s) 
by b(s). From equation (|2.4|) and the invertibility of A(s), we can conclude that x(s) 
is analytic in a region containing [—1,1]. 

Equation (|2.4p reveals the underlying structure of the solution as a function of s. 
If A(s) and b(s) depend polynomially on s, then (]2 .4[) tells us that x(s) is a rational 
function. Note also that this structure is independent of the particular weight function 
w(s). 

3. Spectral Methods. In this section, we derive the spectral methods we use 
to approximate the solution x(s). We begin with a brief review of the relevant theory 
of orthogonal polynomials, Gaussian quadrature, and Fourier series. We include this 
section primarily for the sake of notation and refer the reader to a standard text 
on orthogonal polynomials [26 for further theoretical details and [15] for a modern 
perspective on computation. 

3.1. Orthogonal Polynomials and Gaussian Quadrature. Let P be the 

space of real polynomials defined on [—1,1], and let P„ C P be the space of polynomials 
of degree at most n. For any p, q in P, we define the inner product as 

(pq) = / p(s)q(s)w(s) ds. (3.1) 



We define a norm on P as \\p\\ L 2 — \J (p 2 ), which is the standard 1? norm for the 
given weight w(s). Let {^k(s)} be the set of polynomials that are orthonormal with 
respect to w(s), i.e. (TTiirj) — 5ij. It is known that {irk(s)} satisfy the three-term 
recurrence relation 

/3& + i7Tfc + i(s) = (s - ak)-Kk(s) - pkTTk-i(s), k = 0, 1, 2, . . . , (3.2) 

with 7r_i(s) = and ttq(s) = 1. If we consider only the first n equations, then we can 
rewrite (|3.2[) as 

STT k (s) = /3 fc 7Tfe_ 1 (s) + a k TT k (s) + /3fc+l7T fc+1 (s), k = 0, 1 ,...,n- 1 . (3.3) 

Setting 7T n (s) = [ttq(s), 7Ti(s), . . . , 7r„_i(s)] T , we can write this conveniently in matrix 
form as 



sir n (s) = J„7T„(s) + (3 n n n (s)e n 



(3.4) 
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where e„ is a vector of zeros with a one in the last entry, and J„ (known as the Jacobi 
matrix) is a symmetric, tridiagonal matrix defined as 



a (3% 



OL n -2 Pn-l 
Pn-l CKn-1 



(3-5) 



The zeros {Aj} of 7r„(s) are the eigenvalues of 3 n and 7r ra (Ai) are the corresponding 
eigenvectors; this follows directly from (|3.4p . Let Q„ be the orthogonal matrix of 
eigenvectors of J„. Then we write the eigenvalue decomposition of J„ as 



(3-6) 



It is known (c.f. [15) ) that the eigenvalues {A^} are the familiar Gaussian quadrature 
points associated with the weight function w(s). The quadrature weight Vi correspond- 
ing to \ is equal to the square of the first component of the eigenvector associated 
with Ai, i.e. 



(3-7) 



The weights {^} are known to be strictly positive. We will use these facts repeatedly 
in the sequel. For an integrable scalar function f(s), we can approximate its integral by 
an n-point Gaussian quadrature rule, which is a weighted sum of function evaluations, 



/i n-l 
f(s)w(s)ds = y2f(\ i )v i +R n (f) 



(3-8) 



If / G P2n-i, then R n {f) = 0; that is to say the degree of exactness of the Gaussian 
quadrature rule is In — 1. We use the notation 



(/)„ = E ^ A >* 



(3-9) 



to denote the Gaussian quadrature rule. This is a discrete approximation to the true 
integral. 

3.2. Fourier Series. The polynomials {^(s)} form an orthonormal basis for 
the Hilbert space 



L 1 = Ll([-1, 1]) = {/ : [-1, 1] - R | ||/|| i2 < ex.} . 
Therefore, any f £ L 2 admits a convergent Fourier series 



f(s) = E (-fa*) 7Ffc ( s ) 



(3.10) 



(3.11) 



fc=0 



The coefficients (fi^k) are called the Fourier coefficients. If we truncate the series 
(|3.1ip after n terms, we are left with a polynomial of degree n — 1 that is the best 
approximation polynomial in the L 2 norm. In other words, if we denote 



Pnf{s)=Y J U^)Ms), 



(3.12) 



fc=0 



SPECTRAL METHODS FOR MATRIX EQUATIONS 



5 



then 

\\f-Pnf\y= mf ||/-P|| L2 . (3.13) 

peVn-l 

In fact, the error made by truncating the series is equal to the sum of squares of the 
neglected coefficients, 

oo 

II/-^/I| 2 l 2 = E^) 2 - ( 3 - 14 ) 

k—n 

These properties of the Fourier series motivate the theory and practice of spectral 
methods. 

We have shown that the each element of the solution x(s) of the parameter- 
ized matrix equation is analytic in a region containing the closed interval [—1,1]. 
Therefore it is continuous and bounded on [—1,1], which implies that Xi (s) e L 2 for 
i = 0, . . . , iV — 1. We can thus write the convergent Fourier expansion for each element 
using vector notation as 

oo 

x(s) = ^2 {xK k )n k {s). (3.15) 

fc=0 

Note that we are abusing the bracket notation here, but this will make further manip- 
ulations very convenient. The computational strategy is to choose a truncation level 
n — 1 and estimate the coefficients of the truncated expansion. 

3.3. Spectral Collocation. The term spectral collocation typically refers to 
the technique of constructing a Lagrange interpolating polynomial through the exact 
solution evaluated at the Gaussian quadrature points. Suppose that Aj, i = 0, . . . , n— 1 
are the Gaussian quadrature points for the weight function w(s). We can construct 
an n — 1 degree polynomial interpolant of the solution through these points as 

n-1 

Xc,n(s) = ^2x(\i)£i(s) = X c l„(s). (3.16) 

i=0 

The vector x(Xi) is the solution to the equation A(Xi)x(Xi) = &(A,). The n — 1 degree 
polynomial £j(s) is the standard Lagrange basis polynomial defined as 

The N x n constant matrix X c (the subscript c is for collocation) has one column for 
each x(Xi), and l n (s) is a vector of the Lagrange basis polynomials. 

By construction, the collocation polynomial x c , n interpolates the true solution 
x(s) at the Gaussian quadrature points. We will use this construction to show the 
connection between the pseudospectral method and the Galerkin method. 

3.4. Pseudospectral Methods. Notice that computing the true coefficients of 
the Fourier expansion of x(s) requires the exact solution. The essential idea of the 
pseudospectral method is to approximate the Fourier coefficients of x(s) by a Gaussian 
quadrature rule. In other words, 



Xp,n 



(s) = ^2(xTT k ) n TTk(s) = X p 7T„(s), (3.18) 



i=0 
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where X p is an N x n constant matrix of the approximated Fourier coefficients; the 
subscript p is for pseudospectral. For clarity, we recall 

n-1 

( X7r k) n = x(Ai)%k(^i)vi- (3.19) 

i=0 

where x(Xi) solves A(Xi)x(Xi) = b(Xi). In general, the number of points in the quadra- 
ture rule need not have any relationship to the order of truncation. However, when 
the number of terms in the truncated series is equal to the number of points in the 
quadrature rule, the pseudospectral approximation is equivalent to the collocation 
approximation. This relationship is well-known, but we include following lemma and 
theorem for use in later proofs. 

Lemma 3.1. Let qo be the first row of Q„, and define D qo = diag(qo). The 
matrices X p and X c are related by the equation X p = X c D qo Q^. 

Proof. Write 

Xp(:,fc) = {xn k ) n 



n-1 



3=0 

~ Cl "^lkn(A,)|| 2 ||7 r „(A i )||2 

= X c D qo Q^(:,fc) 

which implies X p = X c D qn Q^ as required. □ 

Theorem 3.2. The n — 1 degree collocation approximation is equal to the n — 1 
degree pseudospectral approximation using an n-point Gaussian quadrature rule, i.e. 



,(*) = x p , n (s). (3.20) 



for all s. 



Proof. Note that the elements of qo are all non-zero, so D 1 exists. Then lemma 



qo 

13.11 implies X c = X p Q n D~ . Using this change of variables, we can write 

x C: „(s) = X c l„(s) = X p Q„D qo 1 l„(s). (3-21) 

Thus it is sufficient to show that Tr n (s) = Q„D_ o 1 l n (s). Since this is just a vector of 
polynomials with degree at most n — 1, we can do this by multiplying each element 
by each orthonormal basis polynomial up to order n — 1 and integrating. Towards 
this end we define = (l„7r^). 

Using the polynomial exactness of the Gaussian quadrature rule, we compute the 
i, j element of 0. 



n-1 
k=0 

7Tj(A 



||7Tn(Ai)||2 ||7T n (Ai)|| 2 

Qn(0,i)Q„(i,i), 
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which implies that © = D qn Q^. Therefore 

= Q„D qo 1 D qo Q^ 



which completes the proof. □ 

Some refer to the pseudospectral method explicitly as an interpolation method [3] . 
See [20] for an insightful interpretation in terms of a discrete projection. Because of 
this property, we will freely interchange the collocation and pseudospectral approxi- 
mations when convenient in the ensuing analysis. 

The work required to compute the pseudospectral approximation is highly depen- 
dent on the parameterized system. In general, we assume that the computation of 
x(Xi) dominates the work; in other words, the cost of computing Gaussian quadra- 
ture formulas is negligible compared to computing the solution to each linear system. 
Then if each x(Xi) costs 0(N 3 ), the pseudospectral approximation with n terms costs 
0{nN 3 ). 

3.5. Spectral Galerkin. The spectral Galerkin method computes a finite di- 
mensional approximation to x(s) such that each element of the equation residual is 
orthogonal to the approximation space. Define 

r(y,s)=A(s)y(s)-b(s). (3.22) 

The finite dimensional approximation space for each component Xi(s) will be the 
space of polynomials of degree at most n — 1. This space is spanned by the first 
n orthonormal polynomials, i.e. span(7r (s), . . . ,7T n _i(s)) = P n -i- We seek an Re- 
valued polynomial x g , n {s) of maximum degree n — 1 such that 

(n{x g , n )n k ) = 0, t = 0,...,JV- 1, fc = 0,...,n-l, (3.23) 

where ri{x g ^ n ) is the ith component of the residual. We can write equations (|3.23[) in 
matrix notation as 

{r{x g<n )-*l) = (3.24) 

or equivalently 

{Ax g ^l) = . (3.25) 

Since each component of x g ^ n (s) is a polynomial of degree at most n— 1, we can write 
its expansion in {irk(s)} as 

n-l 

%g,n(s) = ^^g,k^k(s) = X ff 7T n (s), (3.26) 

fe=0 

where X s is a constant matrix of size N x n; the subscript g is for Galerkin. Then 
equation (|3.25p becomes 



(3.27) 
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Using the vec notation [TBI Section 4.5], we can rewrite (|3.27p as 

(7r„7r£ <g> A) vec(X ff ) = (tt„ <g> b) . (3.28) 

where vec(X ff ) is an Nn x 1 constant vector equal to the columns of X g stacked on 
top of each other. The constant matrix (^-K n TT^ (g> A) has size Nn x Nn and a distinct 
block structure; the i,j block of size N x N is equal to (■K i iTjA). More explicitly, 



{■KnTrZ <g> A) 



(7r 7r ^4) • ■ • (7ro7T n _iA) 
(7r n _i7r A) • ■ • (7r n _i7r n _i/l)^ 



(3.29) 



Similarly, the zth block of the Nnx 1 vector (7r„ ® b) is equal to (bwi), which is exactly 
the ith Fourier coefficient of b(s). 

Since A(s) is bounded and nonsingular for all s G [—1,1], it is straightforward to 
show that Xg 7 n(s) exists and is unique using the classical Galerkin theorems presented 
and summarized in Brenner and Scott [H Chapter 2]. This implies that X g is unique, 
and since b(s) is arbitrary, we conclude that the matrix (7r„7rJ ® A} is nonsingular 
for all finite truncations n. 

The work required to compute the Galerkin approximation depends on how one 
computes the integrals in equation (|3.28j) . If we assume that the cost of forming the 
system is negligible, then the costly part of the computation is solving the system 
(|3.28p . The size of the matrix (lin-ir 1 ^ ® A) is Nn x Nn, so we expect an opera- 
tion count of 0(N 3 n 3 ), in general. However, many applications beget systems with 
sparsity or exploitable structure that can considerably reduce the required work. 

3.6. Summary. We have discussed two classes of spectral methods: (i) the in- 
terpolatory pseudospectral method which approximates the truncated Fourier series 
of x(s) by using a Gaussian quadrature rule to approximate each Fourier coefficient, 
and (ii) the Galerkin projection method which finds an approximation in a finite 
dimensional subspace such that the residual A(s)x giVi (s) — b(s) is orthogonal to the 
approximation space. In general, the n-term pseudospectral approximation requires n 
solutions of the original parameterized matrix equation (|2.2[) evaluated at the Gaus- 
sian quadrature points, while the Galerkin method requires the solution of the coupled 
linear system of equations (|3.28[) that is n times as large as the original parameter- 
ized matrix equation. A rough operation count for the pseudospectral and Galerkin 
approximations is 0(nN 3 ) and 0(n 3 N 3 ), respectively. 

Before discussing asymptotic error estimates, we first derive some interesting and 
useful connections between these two classes of methods. In particular, we can inter- 
pret each method as a set of functions acting on the infinite Jacobi matrix for the 
weight function w(s); the difference between the methods lies in where each truncates 
the infinite system of equations. 

4. Connections Between Pseudospectral and Galerkin. We begin with 
a useful lemma for representing a matrix of Gauss quadrature integrals in terms of 
functions of the Jacobi matrix. 

Lemma 4.1. Let f(s) be a scalar function analytic in a region containing [— 1, 1]. 
Then (fTz n TZ^) n = /(J„). 
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Proof. We examine the i,j element of the n x n matrix /(J„) 

ej /(J n )ej = ef Q„/(A„)Q^e i 
= qf/(An)qj 



ri-l 



= £/(* 



TTi(Afc) ^"(Afc) 



fc=0 
n-l 



k(A fe )|| 2 ||7r(A fe )|| 2 

= £/(A^(A*)*-i(W 

fc=0 

which completes the proof. □ 

Note that Lemma T4.ll generalizes Theorem 3.4 in [T7]- With this in the arsenal, 
we can prove the following theorem relating pseudospectral to Galerkin. 

Theorem 4.2. The pseudospectral solution is equal to an approximation of the 
Galerkin solution where each integral in equation (|3.28p is approximated by an n-point 
Gauss quadrature formula. In other words, X p solves 

(7r„7r£®A) n vec(X p ) = (7r„®b) n . (4.1) 



Proof. Define the N x n matrix B c = [6(Ao) • • • b(X n -i)]- Using the power series 
expansion of A(s) (equation (|2.3[) K we can write the matrix of each collocation solution 

as 

oo 

A(A fe )=^A i A| (4.2) 

i=0 

for A; = 0, . . . ,n— 1. We collect these into one large block-diagonal system by writing 

(fl A ™ ® vec ( X =) = v °c(B c ). (4.3) 

Let I be the N x N identity matrix. Premultiply (|4.3[) by (D qo ® I) , and by commu- 
tativity of diagonal matrices and the mixed product property, it becomes 

(fl A ™ (D i« ® r ) vec ( x c) = (Dqo ® I)vec(B c ). (4.4) 

Premultiplying (|4.4p by (Q„ ® I), properly inserting (Q^ ® I)(Q n ® I) on the left 
hand side, and using the eigenvalue decomposition (|3.6|) . this becomes 

^ jj, ® A^ (Q„ <g> I)(D qo <g> I)vec(X c ) = (Q„ <g> I)(D qo <g> I)vec(B c ). (4.5) 

But note that Lemma |3 . 1 1 implies 

(Q„ ® I)(D qo ® I)vec(X c ) = vec(X p ). (4.6) 
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Using an argument identical to the proof of Lemma |3. 11 we can write 

(Q„ ® I) (D qo <g> I)vec(B c ) = (tt„ <g> b) n (4.7) 
Finally, using Lemma |4~T1 equation (|4. 5|) becomes 

(7r„7r^ (g) A) n vec(Xp) = (tt„ (g) &)„ . (4.8) 

as required. □ 

Theorem l4.2l begets a corollary giving conditions for equivalence between Galerkin 
and pseudospectral approximations. 

Corollary 4.3. If b(s) contains only polynomials of maximum degree and 
A(s) contains only polynomials of maximum degree 1 (i.e. linear functions of s), then 
x g ,n(s) = x Pt n(s) for n > nib for all s G [—1,1]. 

Proof. The parameterized matrix ir n (s)TT n (s) T £g> A(s) has polynomials of degree 
at most In — 1. Thus, by the polynomial exactness of the Gauss quadrature formulas, 

® A) n = ® A) , (tt„ <g> b) n = (tt„ ® b) . (4.9) 

Therefore X ff = X p , and consequently 

%,n(s) = X s 7r n (s) = X p 7T„(s) = x p ,„(s). (4-10) 

as required. □ 

By taking the transpose of equation (|3.27p and following the steps of the proof of 
theorem 14. 2 1 we get another interesting corollary. 

COROLLARY 4.4. First define A(3 n ) to be the Nn x Nn constant matrix with the 
i,j block of size nxn equal to A(i,j)(J n ). Next define b(3 „) to be the Nnxn constant 
matrix with the ith nxn block equal to bi(3 n ). Then the pseudospectral coefficients 
X p satisfy 

A(J„)vec(Xj) = b(3 n )e , (4.11) 

where eo = [1, 0, . . . , 0] T is an n-vector. 

Theorem 14.21 leads to a fascinating connection between the matrix operators in 
the Galerkin and pseudospectral methods, namely that the matrix in the Galerkin 
system is equal to a submatrix of the matrix from a sufficiently larger pseudospectral 
computation. This is the key to understanding the relationship between the Galerkin 
and pseudospectral approximations. In the following lemma, we denote the first r x r 
principal minor of a matrix M by [M] rxr . 

Lemma 4.5. Let A(s) contain only polynomials of degree at most m a , and let 
b(s) contain only polynomials of degree at most nib- Define 



m = m(n) > max 

Then 



2n- 1 



nib 



(4.12) 



(TVnirl ®A)= [(nrnnl ® A) J NnxNn 

{■K n ®b) = [{■Km®b) m ] Nnxl . 
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Proof. The integrands of the matrix (/K n TZ^ <g) A} are polynomials of degree at 
most In + m a — 2. Therefore they can be integrated exactly with a Gauss quadrature 
rule of order m. A similar argument holds for (7r„ ® b). □ 

Combining Lemma 14.51 with corollary 14. 4[ we get the following proposition relat- 
ing the Galerkin coefficients to the Jacobi matrices for A(s) and b(s) that depend 
polynomially on s. 

Proposition 4.6. Let m, m a , and nrn, be defined as in Lemma \4-5\ Define 
[A] n (3 m ) to be the Nn x Nn constant matrix with the i,j block of size n x n equal to 
[A(i, j)(J m )]nxn for i,j = 0, ...,N — 1. Define [b] n (3 m ) to be the Nn x n constant 
matrix with the ith n x n block equal to [bi(3 m )] nxn for i = l,...,iV. Then the 
Galerkin coefficients X g satisfy 

L4]„(J m )vec(Xj) = [6] n (Jm)eo, (4.13) 

where eg = [1, 0, ... , 0] T is an n-vector. 

Notice that Proposition 14.61 provides a way to compute the exact matrix for the 
Galerkin computation without any symbolic manipulation, but beware that m de- 
pends on both n and the largest degree of polynomial in A(s). Written in this form, 
we have no trouble taking m to infinity, and we arrive at the main theorem of this 
section. 

Theorem 4.7. Using the notation of Proposition \4-6\ and corollary \4-4\ the 
coefficients X g of the n-term Galerkin approximation of the solution x(s) to equation 
(|2.2p satisfy the linear system of equations 

L4]„(Joo)vec(Xj) - [6]„(J oo )e , (4.14) 

where eo = [1, 0, . . . , 0] T is an n-vector. 

Proof. Let A^ ma ^(s) be the truncated power series of A(s) up to order m Q , and let 
M" ib )(s) be the truncated power series of b(s) up to order rub- Since A(s) is analytic 
and bounded away from singularity for all s € [— 1 , 1] , there exists an integer M such 
that A^ ma \s) is also bounded away from singularity for all s G [—1, 1] and all m a > M 
(although the bound may be depend on m a ). Assume that m a > M. 

Define m as in equation (|4. 12|) . Then by Proposition ^. 6i the coefficients Xg'™"'" 11 ' 
of the n-term Galerkin approximation to the solution of the truncated system satisfy 

[A("^] n (J m )vec((X(" l «^)) T ) = [fe("^]„(J m )e . (4.15) 

By the definition of m (equation (|4.12p ). equation (|4.15|) holds for all integers greater 
than some minimum value. Therefore, we can take m — > oo without changing the 
solution at all, i.e. 

[^ (m " ) ], 1 (Joo)vec((X( m -" 1 ^) T ) = [6(" l ^]„(J oo )e . (4.16) 
Next we take m a , nib - * oo to get 

[A^ " ]n(Joo) > [^4]n(Joo) 
[6 ( " lb) ]„(Joo) -> M„(Joo) 

which implies 
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as required. □ 

Theorem 14.71 and corollary 14.41 reveal the fundamental difference between the 
Galerkin and pseudospectral approximations. We put them side-by-side for compari- 
son. 

[A]„(J 00 )vec(Xj) = 6„(J oo )e , A(J„)vec(Xj) = 6(J n )e . (4.18) 

The difference lies in where the truncation occurs. For pseudospectral, the infinite 
Jacobi matrix is first truncated, and then the operator is applied. For Galerkin, the 
operator is applied to the infinite Jacobi matrix, and the resulting system is trun- 
cated. The question that remains is whether it matters. As we will see in the error 
estimates in the next section, the interpolating pseudospectral approximation con- 
verges at a rate comparable to the Galerkin approximation, yet requires considerably 
less computational effort. 

5. Error Estimates. Asymptotic error estimates for polynomial approximation 
are well-established in many contexts, and the theory is now considered classical. Our 
goal is to apply the classical theory to relate the rate of geometric convergence to some 
measure of singularity for the solution. We do not seek the tightest bounds in the 
most appropriate norm as in [6], but instead we offer intuition for understanding 
the asymptotic rate of convergence. We also present a residual error estimate that 
may be more useful in practice. We complement the analysis with two representative 
numerical examples. 

To discuss convergence, we need to choose a norm. In the statements and proofs, 
we will use the standard L? and L°° norms generalized to R^-valued functions. 
Definition 5.1. For a function f : K — > l", define the L 2 and L°° norms as 



\\f\\ L * ■= 



N-l r l 

/ fmw(s)ds (5.1) 

i— n J —I 



li"-rwWi SUP l/*( S )U ( 5 ' 2 ) 
0<i<N-l \ _i< s <i 



With these norms, we can state error estimates for both Galerkin and pseudospec- 
tral methods. 

Theorem 5.2 (Galerkin Asymptotic Error Estimate). Let p* be the sum of 
the semi-axes of the greatest ellipse with foci at ±1 in which Xi(s) is analytic for i = 
0, . . . , A— 1. Then for 1 < p < p* , the asymptotic error in the Galerkin approximation 
is 

\\x~x g , n \\ L2 <Cp- n , (5.3) 

where C is a constant independent of n. 

Proof. We begin with the standard error estimate for the Galerkin method [51 
Section 6.4] in the L? norm, 

\\ x - x 9,n\\ L 2 < C\\x - R n x\\ L 2 . (5.4) 

The constant C is independent of n but depends on the extremes of the bounded 
eigenvalues of A(s). Under the consistency hypothesis, the operator R n is a projection 
operator such that 

\\xi — RnXi\\ L2 — ► 0, n— > oo. (5.5) 
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for i = 0, . . . , N — 1. For our purpose, we let R n x be the expansion of x(s) in terms 
of the Chebyshev polynomials, 



R n x(s) 



n-l 

E 

fe=0 



afeT fe (s) 



where T&(s) is the fcth Chebyshev polynomial, and 



a* ■ 



= — / Xl ( S )T k ( S )(l ~ s 2 )- 1 / 2 ds, c k = 



if k = 
otherwise 



(5.6) 



(5.7) 



for i = 0, . . . , N— 1. Since x(s) is continuous for all s G [—1, 1] and iu(s) is normalized, 
we can bound 



||z - R n x\\ L 2 < v^/V ||x — R n x\\ 



(5.8) 



The Chebyshev series converges uniformly for functions that are continuous on [— 1, 1], 
so we can bound 



| a: - R n x\\ L o 



< 



k—n 
oo 

Ei 

k—n 



afcT fe (s) 



(5.9) 
(5.10) 



since —1 < Tk(s) < 1 for all k. To be sure, the quantity |a^| is the component-wise 
absolute value of the constant vector a^, and the norm || • ||oo is the standard infinity 
norm on WL N . 

Using the classical result stated in Section 3], we have 



lim sup \ &k, 

k — >oo 



,1/fe 



1 



o, 



,N- 1 



(5.11) 



where p* is the sum of the semi- axes of the greatest ellipse with foci at ±1 in which 
Xi(s) is analytic. This implies that asymptotically 



O 



m 



0, 



,7V- 1. 



(5.12) 



for pi < p*. We take p = min^ pi, which suffices to prove the estimate (|5.3I) . □ 

Theorem 15 . 2 1 recalls the well-known fact that the convergence of many polynomial 
approximations (e.g. power series, Fourier series) depend on the size of the region 
in the complex plane in which the function is analytic. Thus, the location of the 
singularity nearest the interval [—1, 1] determines the rate at which the approximation 
converges as one includes higher powers in the polynomial approximation. Next we 
derive a similar result for the pseudospectral approximation using the fact that it 
interpolates x(s) at the Gauss points of the weight function w(s). 

Theorem 5.3 (Pseudospectral Asymptotic Error Estimate). Let p* be the sum 
of the semi-axes of the greatest ellipse with foci at ±1 in which Xi(s) is analytic for 
i = 0, . . . , N — 1. Then for 1 < p < p* , the asymptotic error in the pseudospectral 
approximation is 



< c P - n , 



(5.13) 
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where C is a constant independent of n. 

Proof. Recall that x Ci „(s) is the Lagrange interpolant of x(s) at the Gauss points 
of w(s), and let x c ^ n ^(s) be the ith component of x c .„(s). We will use the result 
from Theorem 4.8] that 

J [Xi(s) - x c . n 4s)) 2 w(s) ds < AEl{ Xi ), (5.14) 

where E n {xj) is the error of the best approximation polynomial in the uniform norm. 
We can, again, bound E n (xi) by the error of the Chebyshev expansion (|5.6[) . Using 
Theorem 13 . 2 1 with equation (|5.14p . 

^p.njl^s — ||^ ^c,n||^2 

< 2VN \\x - R n x\\ LOO . 

The remainder of the proof proceeds exactly as the proof of theorem 15.21 □ 

We have shown, using classical approximation theory, that the interpolating pseu- 
dospectral method and the Galerkin method have the same asymptotic rate of geo- 
metric convergence. This rate of convergence depends on the size of the region in 
the complex plane where the functions x(s) are analytic. The structure of the matrix 
equation reveals at least one singularity that occurs when A(s*) is rank-deficient for 
some s* £ R, assuming the right hand side b(s*) does not fortuitously remove it. 
For a general parameterized matrix, this fact may not be useful. However, for many 
parameterized systems in practice, the range of the parameter is dictated by existence 
and/or stability criteria. The value that makes the system singular is often known 
and has some interpretation in terms of the model. In these cases, one may have an 
upper bound on p, which is the sum of the semi-axes of the ellipse of analyticity, and 
this can be used to estimate the geometric rate of convergence a priori. 

We end this section with a residual error estimate - similar to residual error 
estimates for constant matrix equations - that may be more useful in practice than 
the asymptotic results. 

Theorem 5.4. Define the residual r(y, s) as in equation (|3 . 22|) . and let e(y, s) = 
x(s) — y(s) be the K -valued function representing the error in the approximation 
y(s). Then 

Ci \\r(y)\\ L2 < \\e(y)\\ L2 <C 2 \\r(y)\\ L2 (5.15) 

for some constants C\ and Ci, which are independent ofy(s). 

Proof. Since A(s) is non-singular for all s G [—1,1], we can write 

A-^syiy^) = y(s)-A- 1 (s)b(s) = e(y,s) (5.16) 

so that 

\\e(y)\\l* = (e(yf e(y)) 

= (r T (y)A- T A-ir(y)) 

Since A(s) is bounded, so is A^ 1 (s). Therefore, there exist constants and C| that 
depend only on A(s) such that 

CI (r T (y)r(y)) < (e T (y)e(y)) < C 2 * (r T (y)r(y)) . (5.17) 
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Fig. 6.1. The convergence of the spectral methods applied to equation l|6.1|l . The figure on the 
left shows plots the L 2 error as the order of approximation increases, and the figure on the right 
plots the residual error estimate. The stairstep behavior relates to the fact that xo(s) and xi(s) are 
odd functions over [—1, 1], 



Taking the square root yields the desired result. □ 

Theorem 15.41 states that the L 2 norm of the residual behaves like the L 2 norm of 
the error. In many cases, this residual error may be much easier to compute than the 
true L 2 error. However, as in residual error estimates for constant matrix problems, 
the constants in Theorem 15.41 will be large if the bounds on the eigenvalues of A(s) 
are large. We apply these results in the next section with two numerical examples. 

6. Numerical Examples. We examine two simple examples of spectral meth- 
ods applied to parameterized matrix equations. The first is a 2 x 2 symmetric param- 
eterized matrix, and the second comes from a discretized second order ODE. In both 
cases, we relate the convergence of the spectral methods to the size of the region of 
analyticity and verify this relationship numerically. We also compare the behavior of 
the true error to the behavior of the residual error estimate from theorem 15.41 

To keep the computations simple, we use a constant weight function w(s). The 
corresponding orthonormal polynomials are the normalized Legendre polynomials, 
and the Gauss points are the Gauss-Legendre points. 

6.1. A 2 x 2 Parameterized Matrix Equation. Let e > 0, and consider the 
following parameterized matrix equation 



1 + s s~ 




X (s) 




2 


s 1 




x 1 (s)_ 




1 



(6.1) 



For this case, we can easily compute the exact solution, 

, , 2 - s , , 1 + e - 2s . . 

x °( s ) = T~TZ 2> x i\ s ) = T~~ 2- ( 6 - 2 

1 + e — s z 1 + e — s z 

Both of these functions have a poles at s = ±yT + e, so the sum of the semi-axes of 
the ellipse of analyticity is bounded, i.e. p < yjl + e. Notice that the matrix is linear 
in s, and the right hand side has no dependence on s. Thus, corollary 14.31 implies that 
the Galerkin approximation is equal to the pseudospectral approximation for all n; 
there is no need to solve the system (|3.28[) to compute the Galerkin approximation. 
In figure 16.11 we plot both the true L 2 error and the residual error estimate for four 
values of e. The results confirm the analysis. 
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Fig. 6.2. The convergence of the residual error estimate for the Galerkin and pseudo spectral 
approximations applied to the parameterized matrix equation 116. 81 . 



6.2. A Parameterized Second Order ODE. Consider the second order bound- 



ary value problem 



d ( , .du\ 



t e [0, l] 



«(0) 
u(l) 



where, for e > 0, 



a(s, t) = 1 + 4cos(7rs)(i 2 - i), s e [e, 1]. 



The exact solution is 



u(s, t) 



1 



8 cos(7rs) 



ln(l + 4cos(7rs)(i 2 - t)) 



(6.3) 

(6.4) 
(6.5) 



(6.6) 



(6.7) 



The solution u(s,t) has a singularity at s = and t = 1/2. Notice that we have 
adjusted the range of s to be bounded away from by e. We use a standard piecewise 
linear Galerkin finite element method with 512 elements in the t domain to construct 
a stiffness matrix parameterized by s, i.e. 

(K + co8(tts)K 1 )x(s) = b. (6.8) 

Figure 16.21 shows the convergence of the residual error estimate for both Galerkin 
and pseudospectral approximations as n increases. (Despite having the exact solution 
(16. 7p available, we do not present the decay of the L 2 error; it is dominated entirely 
by the discretization error in the t domain.) As e gets closer to zero, the geometric 
convergence rate of the spectral methods degrades considerably. Also, note that each 
element of the parameterized stiffness matrix is an analytic function of s, but figure 
16.21 verifies that the less expensive pseudospectral approximation converges at the 
same rate as the Galerkin approximation. 

7. Summary and Conclusions. We have presented an application of spectral 
methods to parameterized matrix equations. Such parameterized systems arise in 
many applications. The goal of a spectral method is to construct a global polynomial 
approximation of the R^-valued function that satisfies the parameterized system. 
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We derived two basic spectral methods: (i) the interpolatory pseudospectral 
method, which approximates the coefficients of the truncated Fourier series with Gauss 
quadrature formulas, and (ii) the Galerkin method, which finds an approximation in 
a finite dimensional subspace by requiring that the equation residual be orthogonal to 
the approximation space. The primary work involved in the pseudospectral method 
is solving the parameterized system at a finite set of parameter values, whereas the 
Galerkin method requires the solution of a coupled system of equations many times 
larger than the original parameterized system. 

We showed that one can interpret the differences between these two methods 
as a choice of when to truncate an infinite linear system of equations. Employing 
this relationship we derived conditions under which these two approximations are 
equivalent. In this case, there is no reason to solve the large coupled system of 
equations for the Galerkin approximation. 

Using classical techniques, we presented asymptotic error estimates relating the 
decay of the error to the size of the region of analyticity of the solution; we also 
derived a residual error estimate that may be more useful in practice. We verified the 
theoretical developments with two numerical examples: a 2 x 2 matrix equation and 
a finite element discretization of a parameterized second order ODE. 

The popularity of spectral methods for PDEs stems from their infinite (i.e. ge- 
ometric) order of convergence for smooth functions compared to finite difference 
schemes. We have the same advantage in the case of parameterized matrix equa- 
tions, plus the added bonus that there are no boundary conditions to consider. The 
primary concern for these methods is determining the value of the parameter closest 
to the domain that renders the system singular. 
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